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ABSTRACT 

We consider a thin accretion disc warped due to the Bardeen-Petterson effect, pre- 
senting both analytical and numerical solutions for the situation that the two viscosity 
coefficients vary with radius as power law, with the two power law indices not neces- 
sarily equal. The analytical solutions are compared with numerical ones, showing that 
our new analytical solution is more accurate than previous one, which overestimates 
the inclination changing in the outer disc. Our new analytical solution is appropri- 
ate for moderately warped discs, while for extremely misaligned disc, only numerical 
solution is appropriate. 
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1 INTRODUCTION 

Observational evidences are accumulating that accre- 
tion discs around black holes can be warped. Warped 
accretion discs have been directly observ ed by water 
maser observations in NGC4258 (Mivos hi et al. 
iNeufeld fc Maloneyl 



■ [1995I : 

19951: iHermstein et aLl Il996f ) and 



Circinius galaxy ( Greenhill et al.l 20031 ). The lack of cor- 



relation of radio jets in AGNs and the di s c pla ne of host 
galaxy (iKinnev et all I2OO0I : ISchmitt etlii] |2002| ) can also 
be explained by disc warping. Wu et abl (l2008f ) discussed 
the possibility that double-peaked Balmer lines in AGNs be 
emitted by warped disc. Possible evidence for disc warp- 
ing is also found in X-ray binaries, including the misalign- 
ment bet ween jots and orbital planes in GR O J 1655-40 
(jGreene et al. 200]]: lHi"ellming fc Rupenlll995l). and the pre- 
cessing of jets in SS433 (iBlundell fc Bowleij|2004l ). 

Theoretically, warpping can be caused by various 
mechanisms, including tidally i nduced warping by a com- 
panion in a binary system (iTerauern fc BertoutI 1 19931 : 
iLarwood etHI Il99d : [Terguem fc BertoutI Il996l), radiation 
driven or self-inducing w a rping, jMaloney et al l Il996|: 
Maloney fc Begelman' Il997l : iMalonev et al.l Il998l: IPringld 



1996, 1997), mag neti call y driven disc warping, l|Lai 



1999I . 



2003l:lPfeiffer fc Lajbooi) and fr ame dragging driven warp- 
ing (jPardeen fc PettersonI 1 19751 ). Herein we consider the 
sha pe of a disc warped by the la st mechanism. 

iBardeen fc PettersonI (|l975h pointed out that, the com- 
bining effect of Lense-Thirring effect and the viscosity within 
the disc cause the inner part of the disc to be aligned 
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with the central black hole, while the outer part of disc re - 
mains tilted, thus resulting in a warped disc. IPringld (|l992l ) 
derived the dyn a mical equation s of s uch a warped disc. 



iScheuer fc Feileij (|l996l . hereafter ISFgfi) analytically solved 
the equation with a first or der approximation, assum ing con- 
stant viscosity coefficients. iLodato fc Pringlg (|2006h numer- 
ically solved the equations , also a ssuming c onstant v iscosity 
coeffic ients. iMartin etHI l|2007l . hereafter IMPT07| ) gener- 
ahzed ISF96l 's analytical solution to the situation t hat the 
viscos ity coeffici ents v aries as power law, and then, iMartinI 
(|2008l . hereafter IMOSI ) used this solution to fit the maser 
observation of NGC4258's disc. 

We carried on a numerical calculation for a warped disc 
with po wer-law varying and compared the results with 
ImptotI 's analytic al soluti on. The importance of t his wo rk 
lies in such a fact: ImptotI 's analytical solution (andlSF96|'s, 
as well) are based on first order approximation, under the 
assumption of a small inclination angle 6 out ^ 1, while the 
real accretion discs can be strongly misaligned 9out ~ 1, e.g., 
the fitting of NGC4258 shows a strong misaligning. A nu- 
merical calculation is needed to tell exactly how the error 
grows. Our calculation shows a prominent deviation between 
analytical solution and exact solution when the inclination 
angle are large, suggesting that the analytical solutions not 
appropriate for study of NGC4258 or other strongly mis- 
ahgned discs. 

We then proposed another way to extrapolate the small 
6out solution to large 6out situation, and thus find a new an- 
alytical solution. The new solution is also compared with 
numerical calculation and proves to be more accurate for 
large 9out situation. We also generalized the analytical so- 
lutions to the situation that z/i and 1^2 have different power 
index. 
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2 THE BASIC SCENARIO AND EQUATIONS 

We u se the assumptions the same as adopted by iPringlel 
The disc is assumed to be a thin one, consisting 
of concentric (but misahgned) circular gas rings. Each ring 
can be totally described with its surface density E, its an- 
gular velocity n, and its radial velocity Vr- Note that $7 is 
a vector, so that it describes both the speed of the rotation 
n = and the orientation of the ring 1 = Cl/Q,. So, the 
state of the disc can be totally described by the distribution 
of the three quantities with radius -R, E = E(_R), 11 = Q,{R), 
and Vr — Vr{R). Each ring will receive viscous torque from 
neighbouring rings whenever the angular velocity 17 changes 
with radius, ^ 7^ 0. Each ring also receive a Lense-Thirring 
torque from the central black hole whenever it is misaligned 
with the black hole. The dynamical equations under such 
assumption are 



= -^(i^EK)' 



surf 



T, 



= Qf 



Where flpre is the Lense-Thirring precession frequency 



17 rj 



ujp/R 



2G3„ 

c^R3 



(2) 



Lsur/ = ELs is the surface density of angular momentum. 
Lb = R^fl is the specific angular momentum, i.e., the angu- 
lar momentum carried by unit mass. Here we use a dot on 
the head to stand for and the prime symbol "'" to stand 

In this work, we use logarithemic coordinate x = 
ln(_R/i?o), where -Ro is an arbitrarily defined length scale, so 
that all the physical quantities shall be written as functions 
of X. The mass of a ring x ^ x + dx is dm = E ■ 2-KRdR — 
2TiTiadx, where annulus density Ea = R?Ti is the mass on 
unit X interval and unit arc angle. The angular momentum 
of the ring is dL = 2-n'Ladx — 2n'EaT-isdx, where annulus 
angular momentum density La — EaLs is the angular mo- 
mentum on unit x interval and unit arc angle. And we de- 
scribe the radial motion of rings with Vx = ji^r, which is 
the X interval the ring moves in unit time. So the disc can 
be describe with (Ea, La, Vx), as functions of x, and the 
evolution of the disc is described with the evolution of the 
functions with time t. In the following we use a dot on the 
and the prime symbol "'" to stand for 



head to stand for 



With the denotation defined above, the equations can 
be written in a simpler form (nevertheless equivalent to the 
previous form). 



Ea = -(EaK)' 

La = (La) ^ + (La) . + (ta) 

\ / adv V / VIS \ / pre 

— + Tvis + ^pre ^ 

T,„ = Ea(i^iJ^'i + f m') 



(3) 



Note that the 
p a 



here means 



instead of 



and 



3 STEADY STATE SOLUTION FOR 
SLIGHTLY MISALIGNED DISC 

Under Keplerian assumption, the disc can be entirely de- 
picted by a distribution of La. Eliminating redundant vari- 
ables, eqs.(|3} can be rewritten as 

'3 



+ 



iai'j 



+ |3(..i,La)'l 

By l-eq.Q, we get the parallel part of the equation. 



(4) 



(5) 

R^-J ■ K-R^"'^''' 
By eq.Q— lxeq.(l5)l, we get the perpendicular part of the 
equation. 

V2 1 

Yr^' 



'3 

-[r^R 

' V2 1 



La) 1' 
1 



+ {v2 (1') ^La^ + np,.e X La 



1' 



(6) 



When the disc is only slightly misaligned, eqsQ can 
be linearized. Taking the z-axis along the direction of 17p,.e, 
we have 1 = l^^^ + lyG.^ -I- IzG.z ~ 1 



\xy, where X^y = 



lx^x-\-lyG,y, when Ix and ly are small enough for their second- 



order term to be neglected. Then 1' — 



1" = 1. 



/"e^ + ;;'ey, and 1" ■ 1 = (1' 



— e + /' e 
0. Using these 



approximations, the two parts of the angular momentum 
equations becomes 



La 

La^xy 



/3 1 \' / 1 
-[r^RT^La) +(^3.i-^Laj 

+3( 



(7) 



Further using SF96 and MPT07's symbol W = Ix + Uy, 
where i - 

La 

LaW : 



'1, the equation become 



-^-..j,LaW'+{^±LaW' 
+ 3(^1^1-^ La) W' + lQpreLaW 



(8) 



It is not surprising that the first part is all the same with that 
for a planary disc. This means for slightly warped disc we can 
find the solution in two steps. In first step the evolution and 
distribution of La are solved, with the misaligning omitted 
and the disc looked upon as planary. In second step the 
inclination at each radius are found, with La already known. 
This two-step method is much easier than finding the exact 
solution. 

To find a steady state solution, we set the left side of 
eqs.® to zero 



/ 3 1 \ ' / 1 

\^^.-^La) +(^3.1-^La 



'^LaW'+{ 



V2 1 



LaW 



(9) 



^?>{v^^La] W' + iQpreLaW 
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The solution of La are simple 
i^iLa = CoR"^^^ + CiR' 



(10) 



where Co and Ci are constants. Ci is connected with the 
condition at inner boundary, and always become unimpor- 
tant when the concerned region are much larger than inner 
radius. So we discard Ci and get 

,5/2 



VlLa = CoR' 
Substituting the La value back, we get 



: 



+ 



-R 



-1/2 



w 



(11) 



(12) 



If vi and V2 vary with radius as power law v\ — 
vu){R/ RoY^ — !/io exp(/?ia::), V2 = i^20 exp(/32a;), the equa- 
tion of W becomes 



'v2oRq 



exp[(l/2 + /32-/3i)a;]W' 



+iuJp exp[(-l/2 - /3i)a;]W^ = 
Physically we have the boundary conditions 



W 
W 





Wo. 



R^O 

R ^ oo 



(13) 



(14) 



Solving ea. (ll3|l under such boundary conditions, we get 
W = fWoo (15) 
where 



f 



-s"Kn{s) 



and 



and 



r(n) 

1/2 + f32-Pi 
I+P2 



I+P2 



(1-z) 



f2oRo 



exp 



(-i±^.) 



and Kn is the nth order modified Bessel function of the 
second kind. The solution reduces to the MPT07 one (see 
eq.(24) therein) when the two viscosity coefficients vary with 
same index /3i = /32 = P, and further to SF96 solution when 
Pi =132= 0. 

By defining the warping radius as 



Rw — 



l/(l+/32) 



\V20 

the parameter s can be written as 

1 + 32 

2 ,^ ti \ 

s ■ 



(16) 



(17) 



I+P2 

Hereafter we always set Rq = R^, i.e., use the warping ra- 
dius as length unit, thus making the problem scale-free, and 
turning the equation into 

2 / 1 + /32. 



-(l-j)exp(^ — xj 



It is easy to see that the warping radius thus defined is 
where the Lense-Thirring precessing timescale and viscosity 
timescale equals 

1 d3 r»2 

^ -'^ w -^w (19) 

UJp U2{R 




Figure 1. Analytical solutions, ly/Woo against Ix/Wca- The solid 
lines: ft = /3i; the long dash lines: ft = /3i -I- 0.1; the short dash 
lines: ft = /3i — 0.1. For each line style the three lines are for 
/3l =0, 1, 2, respectively, from upside to downside. 



We present here the analytical solutions for several sets 
of Pi and P2- The Pi values are 0, 1, 2, respectively, and 
for each Pi we calculated for /32 = /3i , /32 = /3i -I- 0.1, ft = 
Pi — 0.1. The plane of z axis and 1 at infinite radius, lout, 
is set to be the xz plane, so that lout = (sin Sout, 0, cos Sout) 
and Woo = sm 9 out- For each solution we plot in Fig.l the / 
value in the complex plane, which is equivalent to an ly/Woo 
against Ix/Woo plot. In Fig. 2 we plot the absolute value and 
angle of / (divided by 27r) against radius x. The angle of / 
equals the azimuthal angle ip of 1. The absolute value of / is 
\f\ = W^^ ii#f£:7' and equals for small Bout- So Fig.2 
is also against x and ^ against x plots. We divide ip by 
27r so that the value is now the turns 1 has precessed around 
z axis. The very fast growth of ip in the innermost part of 
disc is not important, because 6 is already very small there, 
meaning the disk is almost aligned will black hole spin. 

In the following we call eq. ([15} "solution A" . 



4 NUMERICAL SOLUTION 

We developed a finite differential code to solve the evolution 
of disc. The state of the disc at each time point is represented 
by the La value upon an uniform grid of x (logarithemic 
grid of R). The time differential of La are then evaluated, 
and then the La value at next time point. We used upstream 
differencing for the advective part in the equation. The code 
is designed with ffexibility to solve various physical problems 
by adjusting the initial condition and boundary condition. 

The code can also be used in finding steady-state so- 
lution. If the boundary condition is fixed to the desired 
setting, and the evolution lasts long enough, in principle 
the disc will always arrive at the wanted steady state so- 
lution. However, the computational cost can be enormous, 
due to the large time scale range involved in the system. 
To ensure the solution reached the steady-state value, the 
time T of disc evolution much be at least several times 
larger than the viscosity time scale T > R'^/vi. On the 
other hand, the maximum time step At to keep the algo- 
rithm numerically stable is determined by the time scale 
for angular momentum viscously diffuse over only one grid, 
2 At < where Ax is the grid size. These two con- 
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Figure 2. Analytical solutions, ^ — - against x (the ascending 
lines) and ^ against x (the descending lines). The (/3i, j32) values 
are given in the legend, x is defined as x = ln(JJ/i?w). 



ditions must hold for the whole calculating region Rin to 
Rout- So the number of time steps needed are determined 

^ Ax-H^f-^'^. As an 

•^4 D in — 4 



0, Rout = 10* 



by(#) /(^ 

\ / max \ 

example, supposing /3i — 

Ax = 0.01, we find T > 10*-^, 2At < 10"^^-^, so that 
the time steps needed are dozens of 10^", absolutely unaf- 
fordable. Our way out of this difficulty is artificially add 
a "speeding up" factor K{R) to the evolutionary equation 
eq.im, changing it to 



(20) 



where jr(La) is the time differential of La given in eq.(|4]). 
This new equation leads to the right steady-state solution 
J-{La) = 0, though its intermediate results (the La values 
found before the disc get steady) is physically meaningless. 
We find K{R) — R^~^^ will make the equation converge 
stably and quickly. 

In this work we set a uniform grid of x from Xin ~ —9.2 
to Xout = 9.2 (corresponding to Ri„ ^ 10~* and Rout ~ 10*, 
the latter large enough to nearly infinity), and the space res- 
olution Ax = 0.01. We used a (i^iLa)' = fi^iLa inner bound- 
ary condition, by adding a "ghost grid" at X — Xin ^^X ^ 
and keep La(a;) = e~^^'^~^''-^'^^lia{xi„), in order to imitate a 
planary disc obeying uiLa oc R^^'^ inside of the inner bound- 
ary. At the outer boundary we set a fixed l,a{x out), with 
an inclination angle to black hole spin axis (set as z axis). 
The plane of z axis and Ija{xout) is set to be xz plane. So 
K^out) = {sin 9 out, 0, COS 9 out), or W(x out) = sin Sout. We use 




Figure 3. Comparing the analytical and numerical solutions in 
the ly/WoD against Ix/Woo plot. The /3 values are /3i = /32 = 3/4. 
From upside to downside, the lines are respectively: 1. Numerical 
solution for Bout = 85°, 2. Solution B for Bout = 85°, 3. Numerical 
solution for Bout = 45° , 4. Solution B for Bout = 45° , 5. Solution A 
for all Bout values, and also all the solutions for s'mBout = 0.01. All 
lines for solution A coincides, because solution A keeps ly/Woo 
and Ix/Woo constant for different Bout- All lines for sin9o„t = 
0.01 almost coincides, showing the error is negligible. In this and 
following two figures, we use solid lines for numerical solution, 
long dash lines for solution B, and short dash lines for solution 
A. 



the "solution B" (explained later) as initial condition to save 
computational cost, though the calculation can converge to 
steady state solution from arbitrary initial condition. 

In Fig.3, Fig. 4 and Fig. 5, the numerical solution are 
shown and compared with solution A. As an example, we 
show the results for /3i = /32 = 3/4, and the inclination an- 
gle at outer boundary to be 9out = arcsin(O.Ol), 30°, 85°, or 
equivalently. Woo = 0.01, 0.5, 0.9962. The numerical solu- 
tion and analytical solution A coincides well when the disc 
is only slightly misaligned |Woo| ^ 1, but when the inclina- 
tion angle is large the two solution deviates strongly. So we 
conclude that solution A is not appropriate for large incli- 
nation angle. In the plot of mass distribution, we use 
because analytical solutions predicts E oc R~'^^ (similar as 
in planary disc). The numerically calculated mass distribu- 
tion differs from analytical solution mainly in the vicinity of 
warping radius, showing a dip there. This is natural because 
the warping there bring forth additional angular momentum 
transfer, so that the gas there falls faster than in the planary 
disc, and thus cause a lower density there. 



5 A NEW ANALYTICAL SOLUTION FOR 
NOT SO SLIGHTLY MISALIGNED DISC 

To find a better analytical solution for more strongly mis- 
aligned disc, we define another measure of misaligning V = 
^(cos ip + i sin ip) , where 9 and (p are the inclination angle and 
azimuthal angle of I, respectively. To the first order approxi- 
mation of 9, W and V equals, W — sin 6'(cos ip+isin ip) ~ V. 
So all the equations for W in sec. 3 also holds for V, hence 
we write 



(21) 
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Figure 4. Comparing the analytical and numerical solutions in 
the g ^ ^ against x (the ascending lines) and ^ against x (the 
descending lines) plots. The /3 values are /3i = /32 = 3/4. From 
upside to downside, the ^ ^ lines are respectively: 1. Numerical 
solution for 6out = 85°, 2. Numerical solution for Bout = 45°, 
3. Solution B for all Bout values, and also all the solutions for 
sin Bout = 0.01. 4. Solution A for Bout = 45°, 5. Solution A for 
Bout = 85° , All lines for solution B coincides, because solution A 
keeps ly/Woc and Ix/Wao constant for different Bout- AH lines for 
sinBout = 0.01 almost coincides, showing the error is negligible. 
For the lower line is the numerical solution for Bout = 85°. 
The upper line is the analytical solutions (solution A and B give 
same ip, unvarying with Bout), and the numerical solution for 
Bout = 45° and sin Bout = 0.01 also coincide with this. 



1.4 - 
1.2 - 




I ' ' ' ' ' ' 1 

-2 2 4 6 8 

Figure 5. Comparing the analytical and numerical solutions in 
the /J'^iS against x plot. The solid line: numerical solutions, for 
s'mBout = 0.01, Bout = 45°, and Bout = 85°, resepectively, from 
upside to downside. The short dash line: analytical solution. 

Hereafter we call this "solution B", and ea. (|15p "solution 
A". The two solutions is equivalent for slight misalignment, 
but behave differently when extrapolated to large inclination 
angle 6out- When 9out varies, solution A keeps sinO/ sinOout 
constant at each R, while solution B keeps 9/6out constant. 
Thus solution A causes too quick a decreasing of 6 at the 
outer disc, while solution B gets rid of this backward. 

Solution B is plotted in Fig. 3 and Fig.4 to compare with 
solution A and numerical solutions. In the W/Woo plot, so- 
lution A keeps unchanged with different 9 out, while solution 
B predicts increasing |W|/Woo with increasing 9out; which 
is closer to the numerical results. In the 9/9out ~ x plot, 
solution B keeps unchanged, while solution A predicts a de- 
creasing 9/9out with increasing 9 out, which contradicts the 



numerical results. However, when 9out are so large as 85°, 
even solution B become very inaccurate. On the other hand, 
for very small 9 out, the two analytical solutions are equiva- 
lent and both very accurate. 



6 CONCLUSIONS 

We generalized MPT07's analytical solution of warped ac- 
cretion discs to the situation that the power law index of the 
two viscosity coefficients is not necessarily equal (solution 

A) . We then proposed a new analytical solution (solution 

B) , which is supposed to be more accurate then solution A. 
We also presented the numerical solutions of the dynami- 
cal equations for warped disc. Our comparison between the 
two analytical solutions and the numerical results show that 
solution B is indeed better and is recommendable for mod- 
erately or slightly misaligned disc. For extremely misaligned 
disc, only numerical solution is appropriate. As for the situa- 
tion in NGC4258, MOS's fitting suggested a large inclination 
angle, so that numerical solution is needed for more accurate 
fitting. 
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